!---------------------------------------------------------------
      subroutine output
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use corgan_com_M
      use cindex_com_M
      use numpar_com_M
      use cprplt_com_M, ONLY:
      use cophys_com_M
      use ctemp_com_M, ONLY: fmf, fef, fvf, fbxf, fbyf, fbzf, twx
      use blcom_com_M
      use cplot_com_M
      use netcdf_M
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
!
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer , dimension(8) :: iv
      integer , dimension(4) :: lioinput, lioprint
      integer , dimension(nreg) :: icount
      integer :: jdx, jdy, jdz, ircode, is
      real(double), dimension(1) :: ixi1, ixi2, ieta1, ieta2
      real(double) :: ixi4, ieta4, ixi5, ieta5
      real(double), dimension(8) :: wght
      real(double), dimension(100) :: fpxf, fpyf, fpzf, fpxft, fpyft, fpzft, &
         ucm, vcm, wcm, cm
      real(double), dimension(3,20) :: wsin, wcos
      real(double), dimension(3,12,20) :: tsi, rot
      real(double) :: jxm0, jxm1, jxm2, jxm3, jxm4, jym0, jym1, jym2, jym3, &
         jym4, jzm0, jzm1, jzm2, jzm3, jzm4, jxn0, jxn1, jxn2, jxn3, jxn4, jyn0&
         , jyn1, jyn2, jyn3, jyn4, jzn0, jzn1, jzn2, jzn3, jzn4
      real(double), dimension(20,20) :: immod, remod, immodz, remodz
      real(double) :: aymin, aymax, rfsolar, bxn0, bxn1, bxn2, bxn3, bxn4, bxm0, bxm1, &
         bxm2, bxm3, bxm4, byn0, byn1, byn2, byn3, byn4, bym0, bym1, bym2, bym3&
         , bym4, bzn0, bzn1, bzn2, bzn3, bzn4, bzm0, bzm1, bzm2, bzm3, bzm4, &
         exn0, exn1, exn2, exn3, exn4, exm0, exm1, exm2, exm3, exm4, eyn0, eyn1&
         , eyn2, eyn3, eyn4, eym0, eym1, eym2, eym3, eym4, ezn0, ezn1, ezn2, &
         ezn3, ezn4, ezm0, ezm1, ezm2, ezm3, ezm4
      logical :: expnd, nomore
      character :: name*80
!
      integer :: n
      real(double) :: error, relax, tiny, phibc, dumdt, dumscal, wk
!-----------------------------------------------
!   E x t e r n a l   F u n c t i o n s
!-----------------------------------------------
      integer :: nf_inq_base_pe, nf_set_base_pe, nf_create, nf__create, &
         nf__create_mp, nf_open, nf__open, nf__open_mp, nf_set_fill, nf_redef, &
         nf_enddef, nf__enddef, nf_sync, nf_abort, nf_close, nf_delete, nf_inq&
         , nf_inq_ndims, nf_inq_nvars, nf_inq_natts, nf_inq_unlimdim, &
         nf_def_dim, nf_inq_dimid, nf_inq_dim, nf_inq_dimname, nf_inq_dimlen, &
         nf_rename_dim, nf_inq_att, nf_inq_attid, nf_inq_atttype, nf_inq_attlen&
         , nf_inq_attname, nf_copy_att, nf_rename_att, nf_del_att, &
         nf_put_att_text, nf_get_att_text, nf_put_att_int1, nf_get_att_int1, &
         nf_put_att_int2, nf_get_att_int2, nf_put_att_int, nf_get_att_int, &
         nf_put_att_real, nf_get_att_real, nf_put_att_double, nf_get_att_double&
         , nf_def_var, nf_inq_var, nf_inq_varid, nf_inq_varname, nf_inq_vartype&
         , nf_inq_varndims, nf_inq_vardimid, nf_inq_varnatts, nf_rename_var, &
         nf_copy_var, nf_put_var_text, nf_get_var_text, nf_put_var_int1, &
         nf_get_var_int1, nf_put_var_int2, nf_get_var_int2, nf_put_var_int, &
         nf_get_var_int, nf_put_var_real, nf_get_var_real, nf_put_var_double, &
         nf_get_var_double, nf_put_var1_text, nf_get_var1_text, &
         nf_put_var1_int1, nf_get_var1_int1, nf_put_var1_int2, nf_get_var1_int2&
         , nf_put_var1_int, nf_get_var1_int, nf_put_var1_real, nf_get_var1_real&
         , nf_put_var1_double, nf_get_var1_double, nf_put_vara_text, &
         nf_get_vara_text, nf_put_vara_int1, nf_get_vara_int1, nf_put_vara_int2&
         , nf_get_vara_int2, nf_put_vara_int, nf_get_vara_int, nf_put_vara_real&
         , nf_get_vara_real, nf_put_vara_double, nf_get_vara_double, &
         nf_put_vars_text, nf_get_vars_text, nf_put_vars_int1, nf_get_vars_int1&
         , nf_put_vars_int2, nf_get_vars_int2, nf_put_vars_int, nf_get_vars_int&
         , nf_put_vars_real, nf_get_vars_real, nf_put_vars_double, &
         nf_get_vars_double, nf_put_varm_text, nf_get_varm_text, &
         nf_put_varm_int1, nf_get_varm_int1, nf_put_varm_int2, nf_get_varm_int2&
         , nf_put_varm_int, nf_get_varm_int, nf_put_varm_real, nf_get_varm_real&
         , nf_put_varm_double, nf_get_varm_double, nccre, ncopn, ncddef, ncdid&
         , ncvdef, ncvid, nctlen, ncsfil
      logical :: nf_issyserr
      character :: nf_inq_libvers*80, nf_strerror*80
!-----------------------------------------------
!
!	LEGENDA
!
!	output switches: 1=on, 0=off
!
!	2 qdnc
!	3 num
!	4 pot
!	5 currents
!	6 B fileds
!	7 E fields
!	8 species quantities (density, velocity, pressure tensor)
!	10 1D fourier analysis
!	11 2D fourier analysis
!	20 particle data
!     compute Ay
!


      call aysub (ay, aymin, aymax, bxn, bzn, iwid, jwid, kwid, ibp1, jbp1, &
         kbp1, dx, dy, dz)
!     compute distance between o points in 2 magnetic islands coalescence
      call opoints(x,bzn,iwid,jwid,kwid,ibp1,jbp1,kbp1,rfsolar)

!
     if (iout(10) /= 0) then
!
!     Calculate fft modes in y
!
      jdx = ibp2
      jdy = jbp2
      jdz = kbp2
!
!     Bx
!
      call modefft1d (ncyc, ibp2, jbp2, kbp2, bxn, iwid, jwid, kwid, jdx, jdy, &
         jdz, a11, a12, a13, 'vec', bxn0, bxn1, bxn2, bxn3, bxn4, bxm0, bxm1, &
         bxm2, bxm3, bxm4)
!
!     By
!
      call modefft1d (ncyc, ibp2, jbp2, kbp2, byn, iwid, jwid, kwid, jdx, jdy, &
         jdz, a11, a12, a13, 'vec', byn0, byn1, byn2, byn3, byn4, bym0, bym1, &
         bym2, bym3, bym4)
!
!     Bz
!
      call modefft1d (ncyc, ibp2, jbp2, kbp2, bzn, iwid, jwid, kwid, jdx, jdy, &
         jdz, a11, a12, a13, 'vec', bzn0, bzn1, bzn2, bzn3, bzn4, bzm0, bzm1, &
         bzm2, bzm3, bzm4)
!
!     Ex
!
      call modefft1d (ncyc, ibp2, jbp2, kbp2, ex, iwid, jwid, kwid, jdx, jdy, &
         jdz, a11, a12, a13, 'vec', exn0, exn1, exn2, exn3, exn4, exm0, exm1, &
         exm2, exm3, exm4)
!
!     Ey      write(*,*)coef2(1:n)
!	stop
!
      call modefft1d (ncyc, ibp2, jbp2, kbp2, ey, iwid, jwid, kwid, jdx, jdy, &
         jdz, a11, a12, a13, 'vec', eyn0, eyn1, eyn2, eyn3, eyn4, eym0, eym1, &
         eym2, eym3, eym4)
!
!     Ez
!
      call modefft1d (ncyc, ibp2, jbp2, kbp2, ez, iwid, jwid, kwid, jdx, jdy, &
         jdz, a11, a12, a13, 'vec', ezn0, ezn1, ezn2, ezn3, ezn4, ezm0, ezm1, &
         ezm2, ezm3, ezm4)
!
!     J12x
!
      call modefft1d (ncyc, ibp2, jbp2, kbp2, jxtilde, iwid, jwid, kwid, jdx, &
         jdy, jdz, a11, a12, a13, 'vec', jxn0, jxn1, jxn2, jxn3, jxn4, jxm0, &
         jxm1, jxm2, jxm3, jxm4)
!
!     J12y
!
      call modefft1d (ncyc, ibp2, jbp2, kbp2, jytilde, iwid, jwid, kwid, jdx, &
         jdy, jdz, a11, a12, a13, 'vec', jyn0, jyn1, jyn2, jyn3, jyn4, jym0, &
         jym1, jym2, jym3, jym4)
!
!     J12z
!
      call modefft1d (ncyc, ibp2, jbp2, kbp2, jztilde, iwid, jwid, kwid, jdx, &
         jdy, jdz, a11, a12, a13, 'vec', jzn0, jzn1, jzn2, jzn3, jzn4, jzm0, &
         jzm1, jzm2, jzm3, jzm4)
!      write(*,*)'JMODES',jzm0,jzm1,jzm2,jzm3,jzm4
end if

if (iout(11) /= 0) then
!
!    1,1 mode for Bx
!
      call modefft2d (ibp2, jbp2, kbp2, bxn, iwid, jwid, kwid, ibar, jbar, kbar&
         , remod, immod, a11, a12, a13, a21, a22, a23)
!
!    1,1 mode for Bz
!
      call modefft2d (ibp2, jbp2, kbp2, bzn, iwid, jwid, kwid, ibar, jbar, kbar&
         , remodz, immodz, a11, a12, a13, a21, a22, a23)

end if
!
!
!     write history variables every time step:
!     write current time to file
!
      scalar = t
      call ncvpt (cdfidh, idtimeh, ncyc, 1, scalar, ircode)
!     write history variables
      scalar = efnrg(nh)
      call ncvpt (cdfidh, idef, ncyc, 1, scalar, ircode)
      scalar = efnrgx(nh)
      call ncvpt (cdfidh, idefx, ncyc, 1, scalar, ircode)
      scalar = efnrgy(nh)
      call ncvpt (cdfidh, idefy, ncyc, 1, scalar, ircode)
      scalar = efnrgz(nh)
      call ncvpt (cdfidh, idefz, ncyc, 1, scalar, ircode)
      scalar = ebnrg(nh)
      call ncvpt (cdfidh, ideb, ncyc, 1, scalar, ircode)
      scalar = ebnrgx(nh)
      call ncvpt (cdfidh, idebx, ncyc, 1, scalar, ircode)
      scalar = ebnrgy(nh)
      call ncvpt (cdfidh, ideby, ncyc, 1, scalar, ircode)
      scalar = ebnrgz(nh)
      call ncvpt (cdfidh, idebz, ncyc, 1, scalar, ircode)
      scalar = ekenrg(nh)
      call ncvpt (cdfidh, ideke, ncyc, 1, scalar, ircode)
      scalar = ekenrgx(nh)
      call ncvpt (cdfidh, idekex, ncyc, 1, scalar, ircode)
      scalar = ekenrgy(nh)
      call ncvpt (cdfidh, idekey, ncyc, 1, scalar, ircode)
      scalar = ekenrgz(nh)
      call ncvpt (cdfidh, idekez, ncyc, 1, scalar, ircode)
      scalar = ekinrg(nh)
      call ncvpt (cdfidh, ideki, ncyc, 1, scalar, ircode)
      scalar = ekinrgx(nh)
      call ncvpt (cdfidh, idekix, ncyc, 1, scalar, ircode)
      scalar = ekinrgy(nh)
      call ncvpt (cdfidh, idekiy, ncyc, 1, scalar, ircode)
      scalar = ekinrgz(nh)
      call ncvpt (cdfidh, idekiz, ncyc, 1, scalar, ircode)
      scalar = eoenrg(nh)
      call ncvpt (cdfidh, ideoe, ncyc, 1, scalar, ircode)
      scalar = eoinrg(nh)
      call ncvpt (cdfidh, ideoi, ncyc, 1, scalar, ircode)
      scalar = charge(nh)
      call ncvpt (cdfidh, idchrg, ncyc, 1, scalar, ircode)
      scalar = citer_maxwell(nh)
      call ncvpt (cdfidh, idxcm, ncyc, 1, scalar, ircode)
      scalar = citer_poisson(nh)
      call ncvpt (cdfidh, idycm, ncyc, 1, scalar, ircode)
      scalar = costcc(nh)
      call ncvpt (cdfidh, idzcm, ncyc, 1, scalar, ircode)
      scalar = aymin
      call ncvpt (cdfidh, idaymin, ncyc, 1, scalar, ircode)
      scalar = aymax
      call ncvpt (cdfidh, idaymax, ncyc, 1, scalar, ircode)
      scalar = rfsolar
      call ncvpt (cdfidh, idrfsolar, ncyc, 1, scalar, ircode)
!
!     history of modes of B
!
if (iout(10) /= 0) then
      scalar = bxn0
      call ncvpt (cdfidh, ibxn0, ncyc, 1, scalar, ircode)
      scalar = bxn1
      call ncvpt (cdfidh, ibxn1, ncyc, 1, scalar, ircode)
      scalar = bxn2
      call ncvpt (cdfidh, ibxn2, ncyc, 1, scalar, ircode)
      scalar = bxn3
      call ncvpt (cdfidh, ibxn3, ncyc, 1, scalar, ircode)
      scalar = bxn4
      call ncvpt (cdfidh, ibxn4, ncyc, 1, scalar, ircode)
      scalar = byn0
      call ncvpt (cdfidh, ibyn0, ncyc, 1, scalar, ircode)
      scalar = byn1
      call ncvpt (cdfidh, ibyn1, ncyc, 1, scalar, ircode)
      scalar = byn2
      call ncvpt (cdfidh, ibyn2, ncyc, 1, scalar, ircode)
      scalar = byn3
      call ncvpt (cdfidh, ibyn3, ncyc, 1, scalar, ircode)
      scalar = byn4
      call ncvpt (cdfidh, ibyn4, ncyc, 1, scalar, ircode)
      scalar = bzn0
      call ncvpt (cdfidh, ibzn0, ncyc, 1, scalar, ircode)
      scalar = bzn1
      call ncvpt (cdfidh, ibzn1, ncyc, 1, scalar, ircode)
      scalar = bzn2
      call ncvpt (cdfidh, ibzn2, ncyc, 1, scalar, ircode)
      scalar = bzn3
      call ncvpt (cdfidh, ibzn3, ncyc, 1, scalar, ircode)
      scalar = bzn4
      call ncvpt (cdfidh, ibzn4, ncyc, 1, scalar, ircode)
!
!     history of nodes of E
!
      scalar = exn0
      call ncvpt (cdfidh, iexn0, ncyc, 1, scalar, ircode)
      scalar = exn1
      call ncvpt (cdfidh, iexn1, ncyc, 1, scalar, ircode)
      scalar = exn2
      call ncvpt (cdfidh, iexn2, ncyc, 1, scalar, ircode)
      scalar = exn3
      call ncvpt (cdfidh, iexn3, ncyc, 1, scalar, ircode)
      scalar = exn4
      call ncvpt (cdfidh, iexn4, ncyc, 1, scalar, ircode)
      scalar = eyn0
      call ncvpt (cdfidh, ieyn0, ncyc, 1, scalar, ircode)
      scalar = eyn1
      call ncvpt (cdfidh, ieyn1, ncyc, 1, scalar, ircode)
      scalar = eyn2
      call ncvpt (cdfidh, ieyn2, ncyc, 1, scalar, ircode)
      scalar = eyn3
      call ncvpt (cdfidh, ieyn3, ncyc, 1, scalar, ircode)
      scalar = eyn4
      call ncvpt (cdfidh, ieyn4, ncyc, 1, scalar, ircode)
      scalar = ezn0
      call ncvpt (cdfidh, iezn0, ncyc, 1, scalar, ircode)
      scalar = ezn1
      call ncvpt (cdfidh, iezn1, ncyc, 1, scalar, ircode)
      scalar = ezn2
      call ncvpt (cdfidh, iezn2, ncyc, 1, scalar, ircode)
      scalar = ezn3
      call ncvpt (cdfidh, iezn3, ncyc, 1, scalar, ircode)
      scalar = ezn4
      call ncvpt (cdfidh, iezn4, ncyc, 1, scalar, ircode)
!
!     history of nodes of J12
!
      scalar = jxn0
      call ncvpt (cdfidh, ijxn0, ncyc, 1, scalar, ircode)
      scalar = jxn1
      call ncvpt (cdfidh, ijxn1, ncyc, 1, scalar, ircode)
      scalar = jxn2
      call ncvpt (cdfidh, ijxn2, ncyc, 1, scalar, ircode)
      scalar = jxn3
      call ncvpt (cdfidh, ijxn3, ncyc, 1, scalar, ircode)
      scalar = jxn4
      call ncvpt (cdfidh, ijxn4, ncyc, 1, scalar, ircode)
      scalar = jyn0
      call ncvpt (cdfidh, ijyn0, ncyc, 1, scalar, ircode)
      scalar = jyn1
      call ncvpt (cdfidh, ijyn1, ncyc, 1, scalar, ircode)
      scalar = jyn2
      call ncvpt (cdfidh, ijyn2, ncyc, 1, scalar, ircode)
      scalar = jyn3
      call ncvpt (cdfidh, ijyn3, ncyc, 1, scalar, ircode)
      scalar = jyn4
      call ncvpt (cdfidh, ijyn4, ncyc, 1, scalar, ircode)
      scalar = jzn0
      call ncvpt (cdfidh, ijzn0, ncyc, 1, scalar, ircode)
      scalar = jzn1
      call ncvpt (cdfidh, ijzn1, ncyc, 1, scalar, ircode)
      scalar = jzn2
      call ncvpt (cdfidh, ijzn2, ncyc, 1, scalar, ircode)
      scalar = jzn3
      call ncvpt (cdfidh, ijzn3, ncyc, 1, scalar, ircode)
      scalar = jzn4
      call ncvpt (cdfidh, ijzn4, ncyc, 1, scalar, ircode)
!
!     modes in ky
!
!
!     ky modes of B
!
      scalar = bxm0
      call ncvpt (cdfidh, ibxm0, ncyc, 1, scalar, ircode)
      scalar = bxm1
      call ncvpt (cdfidh, ibxm1, ncyc, 1, scalar, ircode)
      scalar = bxm2
      call ncvpt (cdfidh, ibxm2, ncyc, 1, scalar, ircode)
      scalar = bxm3
      call ncvpt (cdfidh, ibxm3, ncyc, 1, scalar, ircode)
      scalar = bxm4
      call ncvpt (cdfidh, ibxm4, ncyc, 1, scalar, ircode)
      scalar = bym0
      call ncvpt (cdfidh, ibym0, ncyc, 1, scalar, ircode)
      scalar = bym1
      call ncvpt (cdfidh, ibym1, ncyc, 1, scalar, ircode)
      scalar = bym2
      call ncvpt (cdfidh, ibym2, ncyc, 1, scalar, ircode)
      scalar = bym3
      call ncvpt (cdfidh, ibym3, ncyc, 1, scalar, ircode)
      scalar = bym4
      call ncvpt (cdfidh, ibym4, ncyc, 1, scalar, ircode)
      scalar = bzm0
      call ncvpt (cdfidh, ibzm0, ncyc, 1, scalar, ircode)
      scalar = bzm1
      call ncvpt (cdfidh, ibzm1, ncyc, 1, scalar, ircode)
      scalar = bzm2
      call ncvpt (cdfidh, ibzm2, ncyc, 1, scalar, ircode)
      scalar = bzm3
      call ncvpt (cdfidh, ibzm3, ncyc, 1, scalar, ircode)
      scalar = bzm4
      call ncvpt (cdfidh, ibzm4, ncyc, 1, scalar, ircode)
!
!     ky modes of E
!
      scalar = exm0
      call ncvpt (cdfidh, iexm0, ncyc, 1, scalar, ircode)
      scalar = exm1
      call ncvpt (cdfidh, iexm1, ncyc, 1, scalar, ircode)
      scalar = exm2
      call ncvpt (cdfidh, iexm2, ncyc, 1, scalar, ircode)
      scalar = exm3
      call ncvpt (cdfidh, iexm3, ncyc, 1, scalar, ircode)
      scalar = exm4
      call ncvpt (cdfidh, iexm4, ncyc, 1, scalar, ircode)
      scalar = eym0
      call ncvpt (cdfidh, ieym0, ncyc, 1, scalar, ircode)
      scalar = eym1
      call ncvpt (cdfidh, ieym1, ncyc, 1, scalar, ircode)
      scalar = eym2
      call ncvpt (cdfidh, ieym2, ncyc, 1, scalar, ircode)
      scalar = eym3
      call ncvpt (cdfidh, ieym3, ncyc, 1, scalar, ircode)
      scalar = eym4
      call ncvpt (cdfidh, ieym4, ncyc, 1, scalar, ircode)
      scalar = ezm0
      call ncvpt (cdfidh, iezm0, ncyc, 1, scalar, ircode)
      scalar = ezm1
      call ncvpt (cdfidh, iezm1, ncyc, 1, scalar, ircode)
      scalar = ezm2
      call ncvpt (cdfidh, iezm2, ncyc, 1, scalar, ircode)
      scalar = ezm3
      call ncvpt (cdfidh, iezm3, ncyc, 1, scalar, ircode)
      scalar = ezm4
      call ncvpt (cdfidh, iezm4, ncyc, 1, scalar, ircode)
!
!     ky modes of J12
!
      scalar = jxm0
      call ncvpt (cdfidh, ijxm0, ncyc, 1, scalar, ircode)
      scalar = jxm1
      call ncvpt (cdfidh, ijxm1, ncyc, 1, scalar, ircode)
      scalar = jxm2
      call ncvpt (cdfidh, ijxm2, ncyc, 1, scalar, ircode)
      scalar = jxm3
      call ncvpt (cdfidh, ijxm3, ncyc, 1, scalar, ircode)
      scalar = jxm4
      call ncvpt (cdfidh, ijxm4, ncyc, 1, scalar, ircode)
      scalar = jym0
      call ncvpt (cdfidh, ijym0, ncyc, 1, scalar, ircode)
      scalar = jym1
      call ncvpt (cdfidh, ijym1, ncyc, 1, scalar, ircode)
      scalar = jym2
      call ncvpt (cdfidh, ijym2, ncyc, 1, scalar, ircode)
      scalar = jym3
      call ncvpt (cdfidh, ijym3, ncyc, 1, scalar, ircode)
      scalar = jym4
      call ncvpt (cdfidh, ijym4, ncyc, 1, scalar, ircode)
      scalar = jzm0
      call ncvpt (cdfidh, ijzm0, ncyc, 1, scalar, ircode)
      scalar = jzm1
      call ncvpt (cdfidh, ijzm1, ncyc, 1, scalar, ircode)
      scalar = jzm2
      call ncvpt (cdfidh, ijzm2, ncyc, 1, scalar, ircode)
      scalar = jzm3
      call ncvpt (cdfidh, ijzm3, ncyc, 1, scalar, ircode)
      scalar = jzm4
      call ncvpt (cdfidh, ijzm4, ncyc, 1, scalar, ircode)

endif

      scalar = chcons1(nh)
      call ncvpt (cdfidh, ichcons1, ncyc, 1, scalar, ircode)
      scalar = chcons2(nh)
      call ncvpt (cdfidh, ichcons2, ncyc, 1, scalar, ircode)

if (iout(11) /= 0) then
      pmod = remod
      j = 21
      i = 21
      call writemat (cdfidh, idrmod, pmod, ncyc, 20, 20)

      pmod = immod
      j = 21
      i = 21
      call writemat (cdfidh, idimod, pmod, ncyc, 20, 20)

      pmod = remodz
      j = 21
      i = 21
      call writemat (cdfidh, idrmodz, pmod, ncyc, 20, 20)

      pmod = immodz
      j = 21
      i = 21
      call writemat (cdfidh, idimodz, pmod, ncyc, 20, 20)
endif
!
!     check if it is time to plot other data:
!
!           if(numit.ge.100) go to 1
      if (ncyc /= 1) then
         if (t + 1.E-10 < tout) return
         ixto = ixto + 1
!           make sure plotting interval is nonzero
!           use the last nonzero plotting interval:
         if (dto(ixto) == 0.) ixto = ixto - 1
         tout = tout + dto(ixto)
      endif
!
      iplot = iplot + 1
!
!     write current time and cycle number to file
!
      scalar = t
      call ncvpt (cdfid, idtime, iplot, 1, scalar, ircode)
      call ncvpt (cdfid, idcyc, iplot, 1, ncyc, ircode)
      call ncvpt (cdfidp, idtimep, iplot, 1, scalar, ircode)
      call ncvpt (cdfidp, idcycp, iplot, 1, ncyc, ircode)

!
!     write other grid variables as required:
!     grid
      do ijk = 1, ibp2*jbp2*kbp2
         plotx(ijk) = x(ijk)
         ploty(ijk) = y(ijk)
         plotz(ijk) = z(ijk)
      end do
!
      call writevar (cdfid, idnx, plotx, iplot, ibp2, jbp2, kbp2)
      call writevar (cdfid, idny, ploty, iplot, ibp2, jbp2, kbp2)
      call writevar (cdfid, idnz, plotz, iplot, ibp2, jbp2, kbp2)
!        needed for particles also
      call writevar (cdfidp, idnxp, plotx, iplot, ibp2, jbp2, kbp2)
      call writevar (cdfidp, idnyp, ploty, iplot, ibp2, jbp2, kbp2)
      call writevar (cdfidp, idnzp, plotz, iplot, ibp2, jbp2, kbp2)
!
!     density
      if (iout(2) /= 0) then
         do ijk = 1, ibp2*jbp2*kbp2
            plotx(ijk) = qdnc(ijk)
         end do
!
         call writevar (cdfid, idrho, plotx, iplot, ibp2, jbp2, kbp2)
      endif
!
!     number
      if (iout(3) /= 0) then
         do is = 1, min(2,nsp)
            do ijk = 1, ibp2*jbp2*kbp2
               plotx(ijk) = number(ijk,is)
            end do
!
            call writevar (cdfid, idnum(is), plotx, iplot, ibp2, jbp2, kbp2)
         end do
      endif
!
!     potential
!
      if (iout(4) /= 0) then
         do ijk = 1, ibp2*jbp2*kbp2
            plotx(ijk) = potavp(ijk)
         end do
!
         call writevar (cdfid, idpot, plotx, iplot, ibp2, jbp2, kbp2)
      endif
!
!     current
      if (iout(5) /= 0) then
         do ijk = 1, ibp2*jbp2*kbp2
            plotx(ijk) = jxtilde(ijk)
            ploty(ijk) = jytilde(ijk)
            plotz(ijk) = jztilde(ijk)
         end do
!
         call writevar (cdfid, idjx, plotx, iplot, ibp2, jbp2, kbp2)
         call writevar (cdfid, idjy, ploty, iplot, ibp2, jbp2, kbp2)
         call writevar (cdfid, idjz, plotz, iplot, ibp2, jbp2, kbp2)
!
         do ijk = 1, ibp2*jbp2*kbp2
            plotx(ijk) = jx0(ijk)/(vvol(ijk)+1.E-10)
            ploty(ijk) = jy0(ijk)/(vvol(ijk)+1.E-10)
            plotz(ijk) = jz0(ijk)/(vvol(ijk)+1.E-10)
         end do
!
         call writevar (cdfid, idj0x, plotx, iplot, ibp2, jbp2, kbp2)
         call writevar (cdfid, idj0y, ploty, iplot, ibp2, jbp2, kbp2)
         call writevar (cdfid, idj0z, plotz, iplot, ibp2, jbp2, kbp2)
      endif
!
!     magnetic field
      if (iout(6) /= 0) then
         do ijk = 1, ibp2*jbp2*kbp2
            plotx(ijk) = bxn(ijk)
            ploty(ijk) = byn(ijk)
            plotz(ijk) = bzn(ijk)
         end do
!
         call writevar (cdfid, idbx, plotx, iplot, ibp2, jbp2, kbp2)
         call writevar (cdfid, idby, ploty, iplot, ibp2, jbp2, kbp2)
         call writevar (cdfid, idbz, plotz, iplot, ibp2, jbp2, kbp2)
!
         do ijk = 1, ibp2*jbp2*kbp2
            plotx(ijk) = ay(ijk)
         end do
!
         call writevar (cdfid, iday, plotx, iplot, ibp2, jbp2, kbp2)
      endif
!
!     electric field
      if (iout(7) /= 0) then
         do ijk = 1, ibp2*jbp2*kbp2
            plotx(ijk) = ex(ijk)
            ploty(ijk) = ey(ijk)
            plotz(ijk) = ez(ijk)
            !plotx(ijk) = exavp(ijk) !+ ex_const
            !ploty(ijk) = eyavp(ijk) !+ ey_const
            !plotz(ijk) = ezavp(ijk) !+ ez_const
         end do
!
         call writevar (cdfid, idex, plotx, iplot, ibp2, jbp2, kbp2)
         call writevar (cdfid, idey, ploty, iplot, ibp2, jbp2, kbp2)
         call writevar (cdfid, idez, plotz, iplot, ibp2, jbp2, kbp2)
      endif
!
!     velocity field
      if (iout(8) /= 0) then
         do ijk = 1, ibp2*jbp2*kbp2
            plotx(ijk) = densp1(ijK)
            ploty(ijk) = densp2(ijk)
         end do
!
         call writevar (cdfid, iddens1, plotx, iplot, ibp2, jbp2, kbp2)
         call writevar (cdfid, iddens2, ploty, iplot, ibp2, jbp2, kbp2)
!
         do ijk = 1, ibp2*jbp2*kbp2
            plotx(ijk) = uxsp1(ijk)/(densp1(ijk)+1d-10)!jzs(ijk,1)
            ploty(ijk) = uysp1(ijk)/(densp1(ijk)+1d-10)!jzs(ijk,1)
            plotz(ijk) = uzsp1(ijK)/(densp1(ijk)+1d-10)!jzs(ijk,1)
         end do
!
         call writevar (cdfid, idux1, plotx, iplot, ibp2, jbp2, kbp2)
         call writevar (cdfid, iduy1, ploty, iplot, ibp2, jbp2, kbp2)
         call writevar (cdfid, iduz1, plotz, iplot, ibp2, jbp2, kbp2)
!
         do ijk = 1, ibp2*jbp2*kbp2
            plotx(ijk) = uxsp2(ijk)/(densp2(ijk)+1d-10)!jxs(ijk,2)
            ploty(ijk) = uysp2(ijk)/(densp2(ijk)+1d-10)!jys(ijk,2)
            plotz(ijk) = uzsp2(ijk)/(densp2(ijk)+1d-10)!jzs(ijk,2)
         end do
!
         call writevar (cdfid, idux2, plotx, iplot, ibp2, jbp2, kbp2)
         call writevar (cdfid, iduy2, ploty, iplot, ibp2, jbp2, kbp2)
         call writevar (cdfid, iduz2, plotz, iplot, ibp2, jbp2, kbp2)
!
!
         do ijk = 1, ibp2*jbp2*kbp2
            plotx(ijk) = pixysp2(ijk)
            ploty(ijk) = piyzsp2(ijk)
            plotz(ijk) = pixzsp2(ijk)
         end do
!
         call writevar (cdfid, idpixysp2, plotx, iplot, ibp2, jbp2, kbp2)
         call writevar (cdfid, idpiyzsp2, ploty, iplot, ibp2, jbp2, kbp2)
         call writevar (cdfid, idpixzsp2, plotz, iplot, ibp2, jbp2, kbp2)
!
         do ijk = 1, ibp2*jbp2*kbp2
            plotx(ijk) = pixxsp2(ijk)
            ploty(ijk) = piyysp2(ijk)
            plotz(ijk) = pizzsp2(ijk)
         end do
!
         call writevar (cdfid, idpixxsp2, plotx, iplot, ibp2, jbp2, kbp2)
         call writevar (cdfid, idpiyysp2, ploty, iplot, ibp2, jbp2, kbp2)
         call writevar (cdfid, idpizzsp2, plotz, iplot, ibp2, jbp2, kbp2)

!
         do ijk = 1, ibp2*jbp2*kbp2
            plotx(ijk) = pixysp1(ijk)
            ploty(ijk) = piyzsp1(ijk)
            plotz(ijk) = pixzsp1(ijk)
         end do
!
         call writevar (cdfid, idpixysp1, plotx, iplot, ibp2, jbp2, kbp2)
         call writevar (cdfid, idpiyzsp1, ploty, iplot, ibp2, jbp2, kbp2)
         call writevar (cdfid, idpixzsp1, plotz, iplot, ibp2, jbp2, kbp2)
!
         do ijk = 1, ibp2*jbp2*kbp2
            plotx(ijk) = pixxsp1(ijk)
            ploty(ijk) = piyysp1(ijk)
            plotz(ijk) = pizzsp1(ijk)
         end do
!
         call writevar (cdfid, idpixxsp1, plotx, iplot, ibp2, jbp2, kbp2)
         call writevar (cdfid, idpiyysp1, ploty, iplot, ibp2, jbp2, kbp2)
         call writevar (cdfid, idpizzsp1, plotz, iplot, ibp2, jbp2, kbp2)
!
      endif

!     particles:
      if (iout(20) /= 0) then
!
         if (nrg /= 0) then

!
            write (6, *) 'calling writepar'
            call writeint (cdfidp, idnumtot, numtot(1), iplot, nrg)

            call repack (nrg, icount, numtot, ncells, ijkcell, ijkctmp, link, &
               iphead, iphd2, ico, icotmp, px, pxtmp)
            call writepar (cdfidp, idxp, pxtmp, iplot, nptotl)

            call repack (nrg, icount, numtot, ncells, ijkcell, ijkctmp, link, &
               iphead, iphd2, ico, icotmp, py, pxtmp)
            call writepar (cdfidp, idyp, pxtmp, iplot, nptotl)

            call repack (nrg, icount, numtot, ncells, ijkcell, ijkctmp, link, &
               iphead, iphd2, ico, icotmp, pz, pxtmp)
            call writepar (cdfidp, idzp, pxtmp, iplot, nptotl)

            call repack (nrg, icount, numtot, ncells, ijkcell, ijkctmp, link, &
               iphead, iphd2, ico, icotmp, up, pxtmp)
            call writepar (cdfidp, idup, pxtmp, iplot, nptotl)

            call repack (nrg, icount, numtot, ncells, ijkcell, ijkctmp, link, &
               iphead, iphd2, ico, icotmp, vp, pxtmp)
            call writepar (cdfidp, idvp, pxtmp, iplot, nptotl)

            call repack (nrg, icount, numtot, ncells, ijkcell, ijkctmp, link, &
               iphead, iphd2, ico, icotmp, wp, pxtmp)
            call writepar (cdfidp, idwp, pxtmp, iplot, nptotl)

            call repack (nrg, icount, numtot, ncells, ijkcell, ijkctmp, link, &
               iphead, iphd2, ico, icotmp, qpar, pxtmp)
            call writepar (cdfidp, idqp, pxtmp, iplot, nptotl)
         endif
      endif
!
!     write restart file
   !   call restartout
!
!     sync files:
      call ncsnc (cdfid, ircode)
      call ncsnc (cdfidh, ircode)
      call ncsnc (cdfidp, ircode)
!
      return
      end subroutine output


!
!---------------------------------------------------------------
      subroutine plotclose
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use corgan_com_M
      use cplot_com_M, ONLY: cdfid, cdfidh, cdfidp
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: ircode
      real(double) :: number, mu, lam, lama, lamb, lamc
      real(double), dimension(1) :: ixi1, ixi2, ieta1, ieta2
      real(double) :: ixi4, ieta4, ixi5, ieta5, mupar, jx0, jy0, jz0, jxtilde, &
         jytilde, jztilde, jxs, jys, jzs, jxtildet, jytildet, jztildet, nux, &
         nuy, nuz, mask, iota
      logical :: drift, testpar, adaptg, fields, explicit, ndrft, expnd, fldivb&
         , nomore
      character :: name*80
!-----------------------------------------------
      call ncclos (cdfid, ircode)
      call ncclos (cdfidh, ircode)
      call ncclos (cdfidp, ircode)
      return
      end subroutine plotclose


!
!---------------------------------------------------------------
      subroutine plotinit(ibp2, jbp2, kbp2, nrg)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use corgan_com_M
      use netcdf_M
      use cplot_com_M
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
!
! Routine to create netcdf file and define variable, dimensions, etc.
!
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer  :: ibp2
      integer  :: jbp2
      integer  :: kbp2
      integer  :: nrg
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer , dimension(4) :: dims
      integer :: ircode, ncr, nxdim, nydim, nzdim, nxdimp, nydimp, nzdimp, &
         npdim, moddim, moddim2, nrgdim, is
      real(double) :: number, mu, lam, lama, lamb, lamc
      real(double), dimension(1) :: ixi1, ixi2, ieta1, ieta2
      real(double) :: ixi4, ieta4, ixi5, ieta5, mupar, jx0, jy0, jz0, jxtilde, &
         jytilde, jztilde, jxs, jys, jzs, jxtildet, jytildet, jztildet, nux, &
         nuy, nuz, mask, iota
      logical :: drift, testpar, adaptg, fields, explicit, ndrft, expnd, fldivb&
         , nomore
      character, dimension(nsp) :: label*7
!-----------------------------------------------
!   E x t e r n a l   F u n c t i o n s
!-----------------------------------------------
      integer :: nf_inq_base_pe, nf_set_base_pe, nf_create, nf__create, &
         nf__create_mp, nf_open, nf__open, nf__open_mp, nf_set_fill, nf_redef, &
         nf_enddef, nf__enddef, nf_sync, nf_abort, nf_close, nf_delete, nf_inq&
         , nf_inq_ndims, nf_inq_nvars, nf_inq_natts, nf_inq_unlimdim, &
         nf_def_dim, nf_inq_dimid, nf_inq_dim, nf_inq_dimname, nf_inq_dimlen, &
         nf_rename_dim, nf_inq_att, nf_inq_attid, nf_inq_atttype, nf_inq_attlen&
         , nf_inq_attname, nf_copy_att, nf_rename_att, nf_del_att, &
         nf_put_att_text, nf_get_att_text, nf_put_att_int1, nf_get_att_int1, &
         nf_put_att_int2, nf_get_att_int2, nf_put_att_int, nf_get_att_int, &
         nf_put_att_real, nf_get_att_real, nf_put_att_double, nf_get_att_double&
         , nf_def_var, nf_inq_var, nf_inq_varid, nf_inq_varname, nf_inq_vartype&
         , nf_inq_varndims, nf_inq_vardimid, nf_inq_varnatts, nf_rename_var, &
         nf_copy_var, nf_put_var_text, nf_get_var_text, nf_put_var_int1, &
         nf_get_var_int1, nf_put_var_int2, nf_get_var_int2, nf_put_var_int, &
         nf_get_var_int, nf_put_var_real, nf_get_var_real, nf_put_var_double, &
         nf_get_var_double, nf_put_var1_text, nf_get_var1_text, &
         nf_put_var1_int1, nf_get_var1_int1, nf_put_var1_int2, nf_get_var1_int2&
         , nf_put_var1_int, nf_get_var1_int, nf_put_var1_real, nf_get_var1_real&
         , nf_put_var1_double, nf_get_var1_double, nf_put_vara_text, &
         nf_get_vara_text, nf_put_vara_int1, nf_get_vara_int1, nf_put_vara_int2&
         , nf_get_vara_int2, nf_put_vara_int, nf_get_vara_int, nf_put_vara_real&
         , nf_get_vara_real, nf_put_vara_double, nf_get_vara_double, &
         nf_put_vars_text, nf_get_vars_text, nf_put_vars_int1, nf_get_vars_int1&
         , nf_put_vars_int2, nf_get_vars_int2, nf_put_vars_int, nf_get_vars_int&
         , nf_put_vars_real, nf_get_vars_real, nf_put_vars_double, &
         nf_get_vars_double, nf_put_varm_text, nf_get_varm_text, &
         nf_put_varm_int1, nf_get_varm_int1, nf_put_varm_int2, nf_get_varm_int2&
         , nf_put_varm_int, nf_get_varm_int, nf_put_varm_real, nf_get_varm_real&
         , nf_put_varm_double, nf_get_varm_double, nccre, ncopn, ncddef, ncdid&
         , ncvdef, ncvid, nctlen, ncsfil
      logical :: nf_issyserr
      character :: nf_inq_libvers*80, nf_strerror*80
!-----------------------------------------------
!
!     initialize counter for plots
      iplot = 0
!
!     create netcdf files
      cdfid = nccre('Cel.dat',ncclob,ircode)
      cdfidh = nccre('Celh.dat',ncclob,ircode)
      cdfidp = nccre('Celp.dat',ncclob,ircode)
!
! put name from input file in as global attribute
      ncr = 80
      call ncaptc (cdfid, ncglobal, 'name', ncchar, 80, nome_run, ircode)
      call ncaptc (cdfidh, ncglobal, 'name', ncchar, 80, nome_run, ircode)
      call ncaptc (cdfidp, ncglobal, 'name', ncchar, 80, nome_run, ircode)
! put date, time as global attributes
      call ncaptc (cdfid, ncglobal, 'date', ncchar, 10, ddate, ircode)
      call ncaptc (cdfid, ncglobal, 'hour', ncchar, 10, hour, ircode)
      call ncaptc (cdfidh, ncglobal, 'date', ncchar, 10, ddate, ircode)
      call ncaptc (cdfidh, ncglobal, 'hour', ncchar, 10, hour, ircode)
      call ncaptc (cdfidp, ncglobal, 'date', ncchar, 10, ddate, ircode)
      call ncaptc (cdfidp, ncglobal, 'hour', ncchar, 10, hour, ircode)
!
! define dimensions & get netcdf dimension ids (unlimited must be first?)
      timdimh = ncddef(cdfidh,'time',ncunlim,ircode)
      timdim = ncddef(cdfid,'time',ncunlim,ircode)
      timdimp = ncddef(cdfidp,'time',ncunlim,ircode)
      nxdim = ncddef(cdfid,'nx',ibp2,ircode)
      nydim = ncddef(cdfid,'ny',jbp2,ircode)
      nzdim = ncddef(cdfid,'nz',kbp2,ircode)
      nxdimp = ncddef(cdfidp,'nx',ibp2,ircode)
      nydimp = ncddef(cdfidp,'ny',jbp2,ircode)
      nzdimp = ncddef(cdfidp,'nz',kbp2,ircode)
      npdim = ncddef(cdfidp,'np',npart,ircode)
      moddim = ncddef(cdfidh,'mx',20,ircode)
      moddim2 = ncddef(cdfidh,'my',20,ircode)
      if (nrg /= 0) nrgdim = ncddef(cdfidp,'nrg',nrg,ircode)
!
! define time variable & get id
      dims(1) = timdim
      idtime = ncvdef(cdfid,'t',ncfloat,1,dims,ircode)
      call ncaptc (cdfid, idtime, 'data_type', ncchar, 4, 'hide', ircode)
      dims(1) = timdimh
      idtimeh = ncvdef(cdfidh,'t',ncfloat,1,dims,ircode)
      dims(1) = timdimp
      idtimep = ncvdef(cdfidp,'t',ncfloat,1,dims,ircode)
!
! define cyle variable and get id
      dims(1) = timdim
      idcyc = ncvdef(cdfid,'ncyc',nclong,1,dims,ircode)
      call ncaptc (cdfid, idcyc, 'data_type', ncchar, 4, 'hide', ircode)
      dims(1) = timdim
      idcycp = ncvdef(cdfidp,'ncyc',nclong,1,dims,ircode)
      call ncaptc (cdfidp, idcycp, 'data_type', ncchar, 4, 'hide', ircode)
!
! define history variables
      dims(1) = timdimh
      idef = ncvdef(cdfidh,'efnrg',ncfloat,1,dims,ircode)
      idefx = ncvdef(cdfidh,'efnrgx',ncfloat,1,dims,ircode)
      idefy = ncvdef(cdfidh,'efnrgy',ncfloat,1,dims,ircode)
      idefz = ncvdef(cdfidh,'efnrgz',ncfloat,1,dims,ircode)
      ideb = ncvdef(cdfidh,'ebnrg',ncfloat,1,dims,ircode)
      idebx = ncvdef(cdfidh,'ebnrgx',ncfloat,1,dims,ircode)
      ideby = ncvdef(cdfidh,'ebnrgy',ncfloat,1,dims,ircode)
      idebz = ncvdef(cdfidh,'ebnrgz',ncfloat,1,dims,ircode)
      ideke = ncvdef(cdfidh,'ekenrg',ncfloat,1,dims,ircode)
      idekex = ncvdef(cdfidh,'ekenrgx',ncfloat,1,dims,ircode)
      idekey = ncvdef(cdfidh,'ekenrgy',ncfloat,1,dims,ircode)
      idekez = ncvdef(cdfidh,'ekenrgz',ncfloat,1,dims,ircode)
      ideki = ncvdef(cdfidh,'ekinrg',ncfloat,1,dims,ircode)
      idekix = ncvdef(cdfidh,'ekinrgx',ncfloat,1,dims,ircode)
      idekiy = ncvdef(cdfidh,'ekinrgy',ncfloat,1,dims,ircode)
      idekiz = ncvdef(cdfidh,'ekinrgz',ncfloat,1,dims,ircode)
      ideoe = ncvdef(cdfidh,'eoenrg',ncfloat,1,dims,ircode)
      ideoi = ncvdef(cdfidh,'eoinrg',ncfloat,1,dims,ircode)
      idchrg = ncvdef(cdfidh,'charge',ncfloat,1,dims,ircode)
      idxcm = ncvdef(cdfidh,'itmaxw',ncfloat,1,dims,ircode)
      idycm = ncvdef(cdfidh,'itpois',ncfloat,1,dims,ircode)
      idzcm = ncvdef(cdfidh,'costcc',ncfloat,1,dims,ircode)
      idaymin = ncvdef(cdfidh,'aymin',ncfloat,1,dims,ircode)
      idaymax = ncvdef(cdfidh,'aymax',ncfloat,1,dims,ircode)
      idrfsolar = ncvdef(cdfidh,'rfsolar',ncfloat,1,dims,ircode)
if (iout(10) /= 0) then
      ibxn0 = ncvdef(cdfidh,'bxn0',ncfloat,1,dims,ircode)
      ibxn1 = ncvdef(cdfidh,'bxn1',ncfloat,1,dims,ircode)
      ibxn2 = ncvdef(cdfidh,'bxn2',ncfloat,1,dims,ircode)
      ibxn3 = ncvdef(cdfidh,'bxn3',ncfloat,1,dims,ircode)
      ibxn4 = ncvdef(cdfidh,'bxn4',ncfloat,1,dims,ircode)
      ibyn0 = ncvdef(cdfidh,'byn0',ncfloat,1,dims,ircode)
      ibyn1 = ncvdef(cdfidh,'byn1',ncfloat,1,dims,ircode)
      ibyn2 = ncvdef(cdfidh,'byn2',ncfloat,1,dims,ircode)
      ibyn3 = ncvdef(cdfidh,'byn3',ncfloat,1,dims,ircode)
      ibyn4 = ncvdef(cdfidh,'byn4',ncfloat,1,dims,ircode)
      ibzn0 = ncvdef(cdfidh,'bzn0',ncfloat,1,dims,ircode)
      ibzn1 = ncvdef(cdfidh,'bzn1',ncfloat,1,dims,ircode)
      ibzn2 = ncvdef(cdfidh,'bzn2',ncfloat,1,dims,ircode)
      ibzn3 = ncvdef(cdfidh,'bzn3',ncfloat,1,dims,ircode)
      ibzn4 = ncvdef(cdfidh,'bzn4',ncfloat,1,dims,ircode)

      iexn0 = ncvdef(cdfidh,'exn0',ncfloat,1,dims,ircode)
      iexn1 = ncvdef(cdfidh,'exn1',ncfloat,1,dims,ircode)
      iexn2 = ncvdef(cdfidh,'exn2',ncfloat,1,dims,ircode)
      iexn3 = ncvdef(cdfidh,'exn3',ncfloat,1,dims,ircode)
      iexn4 = ncvdef(cdfidh,'exn4',ncfloat,1,dims,ircode)
      ieyn0 = ncvdef(cdfidh,'eyn0',ncfloat,1,dims,ircode)
      ieyn1 = ncvdef(cdfidh,'eyn1',ncfloat,1,dims,ircode)
      ieyn2 = ncvdef(cdfidh,'eyn2',ncfloat,1,dims,ircode)
      ieyn3 = ncvdef(cdfidh,'eyn3',ncfloat,1,dims,ircode)
      ieyn4 = ncvdef(cdfidh,'eyn4',ncfloat,1,dims,ircode)
      iezn0 = ncvdef(cdfidh,'ezn0',ncfloat,1,dims,ircode)
      iezn1 = ncvdef(cdfidh,'ezn1',ncfloat,1,dims,ircode)
      iezn2 = ncvdef(cdfidh,'ezn2',ncfloat,1,dims,ircode)
      iezn3 = ncvdef(cdfidh,'ezn3',ncfloat,1,dims,ircode)
      iezn4 = ncvdef(cdfidh,'ezn4',ncfloat,1,dims,ircode)

      ijxn0 = ncvdef(cdfidh,'jxn0',ncfloat,1,dims,ircode)
      ijxn1 = ncvdef(cdfidh,'jxn1',ncfloat,1,dims,ircode)
      ijxn2 = ncvdef(cdfidh,'jxn2',ncfloat,1,dims,ircode)
      ijxn3 = ncvdef(cdfidh,'jxn3',ncfloat,1,dims,ircode)
      ijxn4 = ncvdef(cdfidh,'jxn4',ncfloat,1,dims,ircode)
      ijyn0 = ncvdef(cdfidh,'jyn0',ncfloat,1,dims,ircode)
      ijyn1 = ncvdef(cdfidh,'jyn1',ncfloat,1,dims,ircode)
      ijyn2 = ncvdef(cdfidh,'jyn2',ncfloat,1,dims,ircode)
      ijyn3 = ncvdef(cdfidh,'jyn3',ncfloat,1,dims,ircode)
      ijyn4 = ncvdef(cdfidh,'jyn4',ncfloat,1,dims,ircode)
      ijzn0 = ncvdef(cdfidh,'jzn0',ncfloat,1,dims,ircode)
      ijzn1 = ncvdef(cdfidh,'jzn1',ncfloat,1,dims,ircode)
      ijzn2 = ncvdef(cdfidh,'jzn2',ncfloat,1,dims,ircode)
      ijzn3 = ncvdef(cdfidh,'jzn3',ncfloat,1,dims,ircode)
      ijzn4 = ncvdef(cdfidh,'jzn4',ncfloat,1,dims,ircode)

      ibxm0 = ncvdef(cdfidh,'bxm0',ncfloat,1,dims,ircode)
      ibxm1 = ncvdef(cdfidh,'bxm1',ncfloat,1,dims,ircode)
      ibxm2 = ncvdef(cdfidh,'bxm2',ncfloat,1,dims,ircode)
      ibxm3 = ncvdef(cdfidh,'bxm3',ncfloat,1,dims,ircode)
      ibxm4 = ncvdef(cdfidh,'bxm4',ncfloat,1,dims,ircode)
      ibym0 = ncvdef(cdfidh,'bym0',ncfloat,1,dims,ircode)
      ibym1 = ncvdef(cdfidh,'bym1',ncfloat,1,dims,ircode)
      ibym2 = ncvdef(cdfidh,'bym2',ncfloat,1,dims,ircode)
      ibym3 = ncvdef(cdfidh,'bym3',ncfloat,1,dims,ircode)
      ibym4 = ncvdef(cdfidh,'bym4',ncfloat,1,dims,ircode)
      ibzm0 = ncvdef(cdfidh,'bzm0',ncfloat,1,dims,ircode)
      ibzm1 = ncvdef(cdfidh,'bzm1',ncfloat,1,dims,ircode)
      ibzm2 = ncvdef(cdfidh,'bzm2',ncfloat,1,dims,ircode)
      ibzm3 = ncvdef(cdfidh,'bzm3',ncfloat,1,dims,ircode)
      ibzm4 = ncvdef(cdfidh,'bzm4',ncfloat,1,dims,ircode)

      iexm0 = ncvdef(cdfidh,'exm0',ncfloat,1,dims,ircode)
      iexm1 = ncvdef(cdfidh,'exm1',ncfloat,1,dims,ircode)
      iexm2 = ncvdef(cdfidh,'exm2',ncfloat,1,dims,ircode)
      iexm3 = ncvdef(cdfidh,'exm3',ncfloat,1,dims,ircode)
      iexm4 = ncvdef(cdfidh,'exm4',ncfloat,1,dims,ircode)
      ieym0 = ncvdef(cdfidh,'eym0',ncfloat,1,dims,ircode)
      ieym1 = ncvdef(cdfidh,'eym1',ncfloat,1,dims,ircode)
      ieym2 = ncvdef(cdfidh,'eym2',ncfloat,1,dims,ircode)
      ieym3 = ncvdef(cdfidh,'eym3',ncfloat,1,dims,ircode)
      ieym4 = ncvdef(cdfidh,'eym4',ncfloat,1,dims,ircode)
      iezm0 = ncvdef(cdfidh,'ezm0',ncfloat,1,dims,ircode)
      iezm1 = ncvdef(cdfidh,'ezm1',ncfloat,1,dims,ircode)
      iezm2 = ncvdef(cdfidh,'ezm2',ncfloat,1,dims,ircode)
      iezm3 = ncvdef(cdfidh,'ezm3',ncfloat,1,dims,ircode)
      iezm4 = ncvdef(cdfidh,'ezm4',ncfloat,1,dims,ircode)

      ijxm0 = ncvdef(cdfidh,'jxm0',ncfloat,1,dims,ircode)
      ijxm1 = ncvdef(cdfidh,'jxm1',ncfloat,1,dims,ircode)
      ijxm2 = ncvdef(cdfidh,'jxm2',ncfloat,1,dims,ircode)
      ijxm3 = ncvdef(cdfidh,'jxm3',ncfloat,1,dims,ircode)
      ijxm4 = ncvdef(cdfidh,'jxm4',ncfloat,1,dims,ircode)
      ijym0 = ncvdef(cdfidh,'jym0',ncfloat,1,dims,ircode)
      ijym1 = ncvdef(cdfidh,'jym1',ncfloat,1,dims,ircode)
      ijym2 = ncvdef(cdfidh,'jym2',ncfloat,1,dims,ircode)
      ijym3 = ncvdef(cdfidh,'jym3',ncfloat,1,dims,ircode)
      ijym4 = ncvdef(cdfidh,'jym4',ncfloat,1,dims,ircode)
      ijzm0 = ncvdef(cdfidh,'jzm0',ncfloat,1,dims,ircode)
      ijzm1 = ncvdef(cdfidh,'jzm1',ncfloat,1,dims,ircode)
      ijzm2 = ncvdef(cdfidh,'jzm2',ncfloat,1,dims,ircode)
      ijzm3 = ncvdef(cdfidh,'jzm3',ncfloat,1,dims,ircode)
      ijzm4 = ncvdef(cdfidh,'jzm4',ncfloat,1,dims,ircode)
endif

      ichcons1 = ncvdef(cdfidh,'chnorm2',ncfloat,1,dims,ircode)
      ichcons2 = ncvdef(cdfidh,'chnormM',ncfloat,1,dims,ircode)
!
if (iout(11) /= 0) then
      dims(1) = moddim
      dims(2) = moddim2
      dims(3) = timdimh
      idrmod = ncvdef(cdfidh,'rmod',ncfloat,3,dims,ircode)
      idimod = ncvdef(cdfidh,'imod',ncfloat,3,dims,ircode)
!
      dims(1) = moddim
      dims(2) = moddim2
      dims(3) = timdimh
      idrmodz = ncvdef(cdfidh,'rmodz',ncfloat,3,dims,ircode)
      idimodz = ncvdef(cdfidh,'imodz',ncfloat,3,dims,ircode)
endif
!
! define grid variables x,y,z (if(plots) always write these)
      dims(1) = nxdim
      dims(2) = nydim
      dims(3) = nzdim
      dims(4) = timdim
      idnx = ncvdef(cdfid,'x',ncfloat,4,dims,ircode)
      idny = ncvdef(cdfid,'y',ncfloat,4,dims,ircode)
      idnz = ncvdef(cdfid,'z',ncfloat,4,dims,ircode)
      call ncaptc (cdfid, idnx, 'data_type', ncchar, 4, 'hide', ircode)
      call ncaptc (cdfid, idny, 'data_type', ncchar, 4, 'hide', ircode)
      call ncaptc (cdfid, idnz, 'data_type', ncchar, 4, 'hide', ircode)
! define grid for particles
      dims(1) = nxdimp
      dims(2) = nydimp
      dims(3) = nzdimp
      dims(4) = timdimp
      idnxp = ncvdef(cdfidp,'x',ncfloat,4,dims,ircode)
      idnyp = ncvdef(cdfidp,'y',ncfloat,4,dims,ircode)
      idnzp = ncvdef(cdfidp,'z',ncfloat,4,dims,ircode)
      call ncaptc (cdfidp, idnxp, 'data_type', ncchar, 4, 'hide', ircode)
      call ncaptc (cdfidp, idnyp, 'data_type', ncchar, 4, 'hide', ircode)
      call ncaptc (cdfidp, idnzp, 'data_type', ncchar, 4, 'hide', ircode)
!
! define other grid variables, as required:
!
      dims(1) = nxdim
      dims(2) = nydim
      dims(3) = nzdim
      dims(4) = timdim
!
!     density
      if (iout(2) /= 0) then
         idrho = ncvdef(cdfid,'qdnc',ncfloat,4,dims,ircode)
         call ncaptc (cdfid, idrho, 'data_type', ncchar, 4, 'cell', ircode)
      endif
!
!     number
      if (iout(3) /= 0) then
         do is = 1, min(2,nsp)
            label(1) = 'num-1'
            label(2) = 'num-2'
            idnum(is) = ncvdef(cdfid,label(is),ncfloat,4,dims,ircode)
            call ncaptc (cdfid, idnum(is), 'data_type', ncchar, 4, 'cell', &
               ircode)
         end do
      endif
!
!     potential
      if (iout(4) /= 0) then
         idpot = ncvdef(cdfid,'pot',ncfloat,4,dims,ircode)
         call ncaptc (cdfid, idpot, 'data_type', ncchar, 4, 'cell', ircode)
      endif
!
!     current
      if (iout(5) /= 0) then
         idjx = ncvdef(cdfid,'jxhalf',ncfloat,4,dims,ircode)
         idjy = ncvdef(cdfid,'jyhalf',ncfloat,4,dims,ircode)
         idjz = ncvdef(cdfid,'jzhalf',ncfloat,4,dims,ircode)
         call ncaptc (cdfid, idjx, 'data_type', ncchar, 4, 'grid', ircode)
         call ncaptc (cdfid, idjy, 'data_type', ncchar, 4, 'grid', ircode)
         call ncaptc (cdfid, idjz, 'data_type', ncchar, 4, 'grid', ircode)
         idj0x = ncvdef(cdfid,'curlBx',ncfloat,4,dims,ircode)
         idj0y = ncvdef(cdfid,'curlBy',ncfloat,4,dims,ircode)
         idj0z = ncvdef(cdfid,'curlBz',ncfloat,4,dims,ircode)
         call ncaptc (cdfid, idj0x, 'data_type', ncchar, 4, 'grid', ircode)
         call ncaptc (cdfid, idj0y, 'data_type', ncchar, 4, 'grid', ircode)
         call ncaptc (cdfid, idj0z, 'data_type', ncchar, 4, 'grid', ircode)
      endif
!
!     magnetic field
      if (iout(6) /= 0) then
         idbx = ncvdef(cdfid,'bx',ncfloat,4,dims,ircode)
         idby = ncvdef(cdfid,'by',ncfloat,4,dims,ircode)
         idbz = ncvdef(cdfid,'bz',ncfloat,4,dims,ircode)
         iday = ncvdef(cdfid,'ay',ncfloat,4,dims,ircode)
         call ncaptc (cdfid, idbx, 'data_type', ncchar, 4, 'grid', ircode)
         call ncaptc (cdfid, idby, 'data_type', ncchar, 4, 'grid', ircode)
         call ncaptc (cdfid, idbz, 'data_type', ncchar, 4, 'grid', ircode)
         call ncaptc (cdfid, iday, 'data_type', ncchar, 4, 'grid', ircode)
      endif
!
!     electric field
      if (iout(7) /= 0) then
         idex = ncvdef(cdfid,'ex',ncfloat,4,dims,ircode)
         idey = ncvdef(cdfid,'ey',ncfloat,4,dims,ircode)
         idez = ncvdef(cdfid,'ez',ncfloat,4,dims,ircode)
         call ncaptc (cdfid, idex, 'data_type', ncchar, 4, 'grid', ircode)
         call ncaptc (cdfid, idey, 'data_type', ncchar, 4, 'grid', ircode)
         call ncaptc (cdfid, idez, 'data_type', ncchar, 4, 'grid', ircode)
      endif
!
!     velocity field for two species
      if (iout(8) /= 0) then
         iddens1 = ncvdef(cdfid,'densp1',ncfloat,4,dims,ircode)
         iddens2 = ncvdef(cdfid,'densp2',ncfloat,4,dims,ircode)
         call ncaptc (cdfid, iddens1, 'data_type', ncchar, 4, 'grid', ircode)
         call ncaptc (cdfid, iddens2, 'data_type', ncchar, 4, 'grid', ircode)
         idux1 = ncvdef(cdfid,'uxsp1',ncfloat,4,dims,ircode)
         iduy1 = ncvdef(cdfid,'uysp1',ncfloat,4,dims,ircode)
         iduz1 = ncvdef(cdfid,'uzsp1',ncfloat,4,dims,ircode)
         call ncaptc (cdfid, idux1, 'data_type', ncchar, 4, 'grid', ircode)
         call ncaptc (cdfid, iduy1, 'data_type', ncchar, 4, 'grid', ircode)
         call ncaptc (cdfid, iduz1, 'data_type', ncchar, 4, 'grid', ircode)
         idux2 = ncvdef(cdfid,'uxsp2',ncfloat,4,dims,ircode)
         iduy2 = ncvdef(cdfid,'uysp2',ncfloat,4,dims,ircode)
         iduz2 = ncvdef(cdfid,'uzsp2',ncfloat,4,dims,ircode)
         call ncaptc (cdfid, idux2, 'data_type', ncchar, 4, 'grid', ircode)
         call ncaptc (cdfid, iduy2, 'data_type', ncchar, 4, 'grid', ircode)
         call ncaptc (cdfid, iduz2, 'data_type', ncchar, 4, 'grid', ircode)
         idpixysp1 = ncvdef(cdfid,'pixysp1',ncfloat,4,dims,ircode)
         idpiyzsp1 = ncvdef(cdfid,'piyzsp1',ncfloat,4,dims,ircode)
         idpixzsp1 = ncvdef(cdfid,'pixzsp1',ncfloat,4,dims,ircode)
         call ncaptc (cdfid, idpixysp1, 'data_type', ncchar, 4, 'grid', ircode)
         call ncaptc (cdfid, idpiyzsp1, 'data_type', ncchar, 4, 'grid', ircode)
         call ncaptc (cdfid, idpixzsp1, 'data_type', ncchar, 4, 'grid', ircode)
         idpixxsp1 = ncvdef(cdfid,'pixxsp1',ncfloat,4,dims,ircode)
         idpiyysp1 = ncvdef(cdfid,'piyysp1',ncfloat,4,dims,ircode)
         idpizzsp1 = ncvdef(cdfid,'pizzsp1',ncfloat,4,dims,ircode)
         call ncaptc (cdfid, idpixysp1, 'data_type', ncchar, 4, 'grid', ircode)
         call ncaptc (cdfid, idpiyzsp1, 'data_type', ncchar, 4, 'grid', ircode)
         call ncaptc (cdfid, idpixzsp1, 'data_type', ncchar, 4, 'grid', ircode)
         idpixysp2 = ncvdef(cdfid,'pixysp2',ncfloat,4,dims,ircode)
         idpiyzsp2 = ncvdef(cdfid,'piyzsp2',ncfloat,4,dims,ircode)
         idpixzsp2 = ncvdef(cdfid,'pixzsp2',ncfloat,4,dims,ircode)
         call ncaptc (cdfid, idpixysp2, 'data_type', ncchar, 4, 'grid', ircode)
         call ncaptc (cdfid, idpiyzsp2, 'data_type', ncchar, 4, 'grid', ircode)
         call ncaptc (cdfid, idpixzsp2, 'data_type', ncchar, 4, 'grid', ircode)
         idpixxsp2 = ncvdef(cdfid,'pixxsp2',ncfloat,4,dims,ircode)
         idpiyysp2 = ncvdef(cdfid,'piyysp2',ncfloat,4,dims,ircode)
         idpizzsp2 = ncvdef(cdfid,'pizzsp2',ncfloat,4,dims,ircode)
         call ncaptc (cdfid, idpixysp2, 'data_type', ncchar, 4, 'grid', ircode)
         call ncaptc (cdfid, idpiyzsp2, 'data_type', ncchar, 4, 'grid', ircode)
         call ncaptc (cdfid, idpixzsp2, 'data_type', ncchar, 4, 'grid', ircode)
      endif
!
!     particles:
      if (iout(20) /= 0) then
!  define numtot and id:
         if (nrg /= 0) then
            dims(1) = nrgdim
            dims(2) = timdim
            idnumtot = ncvdef(cdfidp,'numtot',nclong,2,dims,ircode)
            call ncaptc (cdfidp, idnumtot, 'data_type', ncchar, 4, 'hide', &
               ircode)
!
!  define particle positions and id:
            dims(1) = npdim
            dims(2) = timdimp
            idxp = ncvdef(cdfidp,'xp',ncfloat,2,dims,ircode)
            idyp = ncvdef(cdfidp,'yp',ncfloat,2,dims,ircode)
            idzp = ncvdef(cdfidp,'zp',ncfloat,2,dims,ircode)
            idup = ncvdef(cdfidp,'up',ncfloat,2,dims,ircode)
            idvp = ncvdef(cdfidp,'vp',ncfloat,2,dims,ircode)
            idwp = ncvdef(cdfidp,'wp',ncfloat,2,dims,ircode)
            idqp = ncvdef(cdfidp,'qp',ncfloat,2,dims,ircode)
            call ncaptc (cdfidp, idxp, 'data_type', ncchar, 4, 'part', ircode)
            call ncaptc (cdfidp, idyp, 'data_type', ncchar, 4, 'part', ircode)
            call ncaptc (cdfidp, idzp, 'data_type', ncchar, 4, 'part', ircode)
            call ncaptc (cdfidp, idup, 'data_type', ncchar, 4, 'part', ircode)
            call ncaptc (cdfidp, idvp, 'data_type', ncchar, 4, 'part', ircode)
            call ncaptc (cdfidp, idwp, 'data_type', ncchar, 4, 'part', ircode)
            call ncaptc (cdfidp, idqp, 'data_type', ncchar, 4, 'part', ircode)
         endif
!
      endif
!
! end definition section
      call ncendf (cdfid, ircode)
      call ncendf (cdfidh, ircode)
      call ncendf (cdfidp, ircode)
!
      return
      end subroutine plotinit


      subroutine writeint(cdfid, idvar, ivar, iplot, ndim)
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer  :: cdfid
      integer  :: idvar
      integer , intent(in) :: iplot
      integer , intent(in) :: ndim
      integer  :: ivar(ndim)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer , dimension(2) :: istart, icount
      integer :: ircode
!-----------------------------------------------
!
!
! generic netcdf write of particle data assuming start at 1
      istart(1) = 1
      istart(2) = iplot
!
      icount(1) = ndim
      icount(2) = 1
!
      call ncvpt (cdfid, idvar, istart, icount, ivar, ircode)
!
      return
      end subroutine writeint


!---------------------------------------------------------------
!---------------------------------------------------------------
      subroutine writepar(cdfid, idvar, var, iplot, npart)
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer  :: cdfid
      integer  :: idvar
      integer , intent(in) :: iplot
      integer , intent(in) :: npart
      real(4)  :: var(npart)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer , dimension(2) :: istart, icount
      integer :: ircode
!-----------------------------------------------
!
!
! generic netcdf write of particle data assuming start at 1
      istart(1) = 1
      istart(2) = iplot
!
      icount(1) = npart
      icount(2) = 1
!
      call ncvpt (cdfid, idvar, istart, icount, var, ircode)
!
      return
      end subroutine writepar


!---------------------------------------------------------------

      subroutine writevar(cdfid, idvar, var, iplot, nxp, nyp, nzp)
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer  :: cdfid
      integer  :: idvar
      integer , intent(in) :: iplot
      integer , intent(in) :: nxp
      integer , intent(in) :: nyp
      integer , intent(in) :: nzp
      real(4)  :: var(nxp,nyp,nzp)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer , dimension(4) :: istart, icount
      integer :: ircode
!-----------------------------------------------
!

! generic netcdf write of grid quantities assuming start at 1
      istart(1) = 1
      istart(2) = 1
      istart(3) = 1
      istart(4) = iplot
!
      icount(1) = nxp
      icount(2) = nyp
      icount(3) = nzp
      icount(4) = 1
!
      write (6, *) 'writing var'
      call ncvpt (cdfid, idvar, istart, icount, var, ircode)
!
      return
      end subroutine writevar


!---------------------------------------------------------------

      subroutine writemat(cdfid, idvar, var, iplot, nxp, nyp)
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer  :: cdfid
      integer  :: idvar
      integer , intent(in) :: iplot
      integer , intent(in) :: nxp
      integer , intent(in) :: nyp
      real(4)  :: var(nxp,nyp)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer , dimension(3) :: istart, icount
      integer :: ircode
!-----------------------------------------------
!

! generic netcdf write of grid quantities assuming start at 1
      istart(1) = 1
      istart(2) = 1
      istart(3) = iplot
!
      icount(1) = nxp
      icount(2) = nyp
      icount(3) = 1
!
      write (6, *) 'writing mat'
      call ncvpt (cdfid, idvar, istart, icount, var, ircode)
!
      return
      end subroutine writemat


!---------------------------------------------------------------
!
      subroutine repack(nrg, icount, numtot, ncells, ijkcell, ijkctmp, link, &
         iphead, iphd2, ico, icotmp, px, pxtmp)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
!
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: nrg
      integer , intent(in) :: ncells
      integer , intent(inout) :: icount(*)
      integer , intent(in) :: numtot(0:*)
      integer , intent(in) :: ijkcell(*)
      integer , intent(inout) :: ijkctmp(*)
      integer , intent(in) :: link(0:*)
      integer , intent(in) :: iphead(*)
      integer , intent(inout) :: iphd2(*)
      integer , intent(in) :: ico(0:*)
      real(4) , intent(out) :: icotmp(*)
      real(4) , intent(out) :: pxtmp(*)
      real(double) , intent(in) :: px(0:*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: j1, j2, kr, newcell, n, ijk, np, iploc, newcell_nxt
!-----------------------------------------------
!
!
!
!      repack particles:
      icount(1) = 1
      do kr = 2, nrg
         icount(kr) = icount(kr-1) + numtot(kr-1)
      end do
!
!
      newcell = 0
      do n = 1, ncells
         ijk = ijkcell(n)
         np = iphead(ijk)
         if (np == 0) cycle
         newcell = newcell + 1
         ijkctmp(newcell) = ijk
         iphd2(ijk) = np
      end do
!
!
   15 continue
      if (newcell == 0) go to 500
!
      do n = 1, newcell
         ijk = ijkctmp(n)
         np = iphd2(ijk)
         kr = ico(np)
         iploc = icount(kr)
         pxtmp(iploc) = px(np)
         icotmp(iploc) = ico(np)
         icount(kr) = icount(kr) + 1
      end do
!
      newcell_nxt = 0
      do n = 1, newcell
         ijk = ijkctmp(n)
         np = link(iphd2(ijk))
         if (np <= 0) cycle
         newcell_nxt = newcell_nxt + 1
         iphd2(ijk) = np
         ijkctmp(newcell_nxt) = ijk
      end do
      newcell = newcell_nxt
!
      go to 15
!
  500 continue
      return
      end subroutine repack


!_________________________________________________
!
      subroutine repack_object(nsp_object, numtot_object, ncells, ijkcell, &
         link_object, iphead_object, itdim, px_object, py_object, pz_object, &
         pxtmp, pytmp, pztmp)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
!
!
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: nsp_object
      integer , intent(in) :: ncells
      integer , intent(in) :: itdim
      integer , intent(inout) :: numtot_object(0:*)
      integer , intent(in) :: ijkcell(*)
      integer , intent(in) :: link_object(0:*)
      integer , intent(in) :: iphead_object(itdim,*)
      real(double) , intent(in) :: px_object(0:*)
      real(double) , intent(in) :: py_object(0:*)
      real(double) , intent(in) :: pz_object(0:*)
      real(double) , intent(out) :: pxtmp(*)
      real(double) , intent(out) :: pytmp(*)
      real(double) , intent(out) :: pztmp(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: is, icount, n, ijk, np
!-----------------------------------------------
!
!
!
!     repack object particles:
!
      numtot_object(1:nsp_object) = 0
!
      icount = 1
!
      do is = 1, nsp_object
         do n = 1, ncells
            ijk = ijkcell(n)
            np = iphead_object(ijk,is)
   10       continue
            if (np == 0) cycle
            pxtmp(icount) = px_object(np)
            pytmp(icount) = py_object(np)
            pztmp(icount) = pz_object(np)
            numtot_object(is) = numtot_object(is) + 1
            icount = icount + 1
            np = link_object(np)
            go to 10
         end do
      end do
!
      return
      end subroutine repack_object
